#ifndef _INTERPOLATION_H_
#define _INTERPOLATION_H_

#include <iostream>
#include <cmath>
#include <vector>
#include <algorithm>
#include <string>
using namespace std;
const static double tao = 1e-8;

/**
 * @brief Factorial
 * 
*/
int factorial(const int &k){
    int p = 1;
    for (int i = 1; i <= p; i++){
        p = p*i;
    }
    return p;
}


/**
 * @class Function
 * @brief An abstract funtion to get value of f(x) and n-th derivative
 * 
*/
class Function{
public:
    //重载运算符()来计算f()
    virtual double operator() (const double &x) const =0;

    // //f()n-th differential
    // virtual double dev(const double &x, const int &n) const{
    //     if(n==0) return (*this)(x);
    //     else if (n >= 1) {
    //         return ((*this).dev(x+tao,n-1)-(*this).dev(x-tao,n-1))/(2*tao);
    //         }
    //     }
    //导函数
    virtual double diff(const double &x) const{
            return ((*this)(x+tao)-(*this)(x-tao))/(2*tao);
        }
};

/**
 * @class Newton_Inter
 * @brief Distinct xi and f(xi),(0<= i <= n)
 * We use n-th polynomials to interpolate
 * 
*/

class Newton_Inter: public Function{
private:
    //结点列x0,...xn
    vector<double> xx;

    //f0,...fn
    vector<double> fx;

    //得到多项式阶数就是(x0,x1,....xn)中的n
    int n;

public:

    //计算得到的系数
    vector<double> coefft;

    Newton_Inter(vector<double> xx,vector<double> fx):
    xx(xx),fx(fx)
    {
        coefft = fx;
        n = xx.size()-1;
        Newton_Inter::Solve_Coefficient();
        //sort(x.begin(),x.end());
    }

    //按照牛顿插值得到（n+1）个系数
    void Solve_Coefficient(){
        for(int i=0;i<=n-1;i++){
            for(int j=0;j<=n-i-1;j++){
                coefft[n-j] = (coefft[n-j]-coefft[n-j-1])/(xx[n-j]-xx[n-j-i-1]);
            }
        }

        // for(int i=0;i<coefft.size();i++){
        //     cout << coefft[i] << " " ;
        // }
    }


    // //牛顿n次插值多项式
    double operator() (const double &x) const{
        double result =coefft[n];
        for(int i=0;i<= n-1;i++){
            result = coefft[n-i-1]+result*(x-xx[n-i-1]);
        }
        return result;
    }
};

/**
 * @class Hermite
 * @brief  xi and mi-th derivative of f(xi),(0<= i <= n)
 * M = sum of mi
 * We use M-th polynomials to interpolate
 * 
*/

class Hermite : public Function{
private:
    //结点列x0,...xn
    vector<double> xx;
    int nx;
    
    //所有f结点的所有导数
    vector<vector<double>> fx_value;
    vector<vector<double>> fx_value_hmt;

public:
    //计算得到的系数
    vector<double> coefft;

    //hermite插值扩展结点
    vector<double> xx_hmt;

    //多项式阶数
    int np;

    //按照各个点的导数得到有重复点的x0,x1,x1,x2,x2
    Hermite(vector<double> xx,vector<vector<double>> fx_value):
    xx(xx),fx_value(fx_value)
    {
        nx = xx.size();
        for(int i=0;i<nx;i++){
            for(int j=0;j< fx_value[i].size();j++){
                coefft.push_back(fx_value[i][0]);
                xx_hmt.push_back(xx[i]);
                fx_value_hmt.push_back(fx_value[i]);
            }
        }
        np = xx_hmt.size()-1;
        Hermite::Solve_Coefficient();
        
    }
    
    //计算多项式所有系数
    void Solve_Coefficient(){
        np = xx_hmt.size()-1;
        for(int i=0;i<=np-1;i++){
            for(int j=0;j<=np-i-1;j++){
                if(xx_hmt[np-j]==xx_hmt[np-j-i-1]){
                    coefft[np-j] = fx_value_hmt[np-j][i+1]/factorial(i+1);
                }
                else {
                    coefft[np-j] = 
                        (coefft[np-j]-coefft[np-j-1])/(xx_hmt[np-j]-xx_hmt[np-j-i-1]);
                }
            }
        }
    }

    //Hermite n次插值多项式
    double operator() (const double &x) const{
        double result =coefft[np];
        for(int i=0;i<= np-1;i++){
            result = coefft[np-i-1]+result*(x-xx_hmt[np-i-1]);
        }
        return result;
    }
};

#endif
